source('/Users/ben/evo-dispersal/LRS/LRS functions.R')

X<-10
Y<-10
longX<-110
gcK<-100

pop<-init.inds(X*Y*gcK, X, Y)
nbrs<-neighbours.init(X, Y)
for (g in 1:500){
	pop<-repro.glob(pop, X, Y, fec.table, gcK, 0.001, plot=F)
	pop<-disperse(pop, nbrs, X)	
}
equil.samp<-sample(pop[,"P"], 1000)
nbrs<-neighbours.init(longX, Y)
for (g in 1:100){
	pop<-repro.glob(pop, longX, Y, fec.table, gcK, 0.001, plot=F)
	pop<-disperse(pop, nbrs, longX)	
}
re.samp<-pop[order(pop[,"X"]),]
re.samp<-re.samp[(length(re.samp[,1])-1000):length(re.samp[,1]), "P"]

save(equil.samp, re.samp, fec.table, file="/Users/ben/Documents/Papers/Beyond Darwin/Figures/FigureData.RData")


png(filename="/Users/ben/Documents/Papers/Beyond Darwin/Figures/LRS Global Histograms.png", width=10, height=20, units="cm", res=400)
par(mfrow=c(2,1))
par(cex.lab=1.5)
par(mar=c(5,5,4,2))
hist(equil.samp, breaks=seq(-0.1, 1, 0.1), ylab="Frequency", main="Trait values at equilibrium", xlab="", ylim=c(0,200))
text(0, 195, "A", cex=1.5)
hist(re.samp, breaks=seq(-0.1, 1, 0.1), ylab="Frequency", main="Trait values on expanding edge", xlab="Dispersal probability")
text(0,342, "B", cex=1.5)
dev.off()